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Abstract: We propose a general modeling framework for marked Poisson processes observed over 
time or space. The modeling approach exploits the connection of the nonhomogeneous Poisson process 
intensity with a density function. Nonparametric Dirichlet process mixtures for this density, combined 
with nonparametric or semiparametric modeling for the mark distribution, yield flexible prior models 
for the marked Poisson process. In particular, we focus on fully nonparametric model formulations 
that build the mark density and intensity function from a joint nonparametric mixture, and provide 
guidelines for straightforward application of these techniques. A key feature of such models is that they 
can yield flexible inference about the conditional distribution for multivariate marks without requiring 
specification of a complicated dependence scheme. We address issues relating to choice of the Dirichlet 
process mixture kernels, and develop methods for prior specification and posterior simulation for full 
inference about functionals of the marked Poisson process. Moreover, we discuss a method for model 
checking that can be used to assess and compare goodness of fit of different model specifications under 
the proposed framework. The methodology is illustrated with simulated and real data sets. 
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1 Introduction 



Marked point process data, occurring on either spatial or temporal domains, is encountered in 
research for biology, ecology, economics, sociology, and numerous other disciplines. When- 
ever interest lies in the intensity of event occurrences as well as the spatial or temporal distri- 
bution of events, the data analysis problem will involve inference for a non-homogeneous point 
process. Moreover, many applications involve marks - a set of random variables associated 
with each point event - such that the data generating mechanism is characterized as a marked 
point process. In marketing, for example, interest may lie in both the location and intensity of 
purchasing behavior as well as consumer choices, and the data may be modeled as a spatial 
point process with purchase events and product choice marks. As another example, in forestry 
interest often lies in estimating the wood-volume characteristics of a plot of land by under- 
standing the distribution and type of tree in a smaller subplot. Hence, the forest can be modeled 
as a spatial point process with tree events marked by trunk size and tree species. 

Non-homogeneous Poisson processes (NHPPs) play a fundamental role in inference for 
data consisting of point event patterns (e.g., |Guttorp[ [1995] |M0ller and Waagepetersen[ |2004[ ), 



and marked NHPPs provide the natural model extension when the point events are accompanied 
by random marks. One reason for the common usage of Poisson processes is their general 
tractability and the simplicity of the associated data likelihood. In particular, for a NHPP, 
PoP(7^, A), defined on the observation window IZ with intensity A(x) for x E TZ, which is a 
non-negative and locally integrable function for aU bounded B CTZ, the following hold true: 

i. For any such ;B, the number of points in S, A^(i3) ~ Po(A(;B)), where A(i3) = J^X{x.)dx. 
is the NHPP cumulative intensity function. 

ii. Given N{B), the point locations within B are i.i.d. with density A(x) / A(x)(ix. 

Here, Po(/i) denotes the Poisson distribution with mean /i. Although IZ can be of arbitrary 
dimension, we concentrate on the common settings of temporal NHPPs with IZ C R^, or 
spatial NHPPs where 7^ C 

This paper develops Bayesian nonparametric mixtures to model the intensity function of 



NHPPs, and will provide a framework for combining this approach with flexible (nonpara- 
metric or semiparametric) modeling for the associated mark distribution. Since we propose 
fully nonparametric mixture modeling for the point process intensity, but within the context of 
Poisson distributions induced by the NHPP assumption, the nature of our modeling approach 
is semiparametric. We are able to take advantage of the above formulation of the NHPP and 
specify the sampling density /(x) = A(x)/A7^ through a Dirichlet process (DP) mixture model, 
where At^ = ^(^) = j-ji A(x)(ix is the total integrated intensity. Crucially, items i and ii above 
imply that the likelihood for a NHPP generated point pattern {xi, . . . , x^v} <ZTZ factorizes as 



N 



P ({x,,}f=i; A(-)) = P ({x,}f=,; A^, /(■)) cx A^ exp(-A^) J] /(x. 



(1) 



such that the NHPP density, /(■), and integrated intensity, K-ji, can be modeled separately. In 
particular, the DP mixture modeling framework for /(■) allows for data-driven inference about 
non-standard intensity shapes and quantification of the associated uncertainty. 



This approach was originally developed by |Kottas and Sans(5| ( |2007| ) in the context of spatial 
NHPPs with emphasis on extreme value analysis problems, and has also been applied to anal- 



ysis of immunological studies (Ji et al. , 2009) and neuronal data analysis (Kottas and Behseta 



2010[ ). Here, we generalize the mixture model to alternative kernel choices that provide for 
conditionally conjugate models and, in the context of temporal NHPPs, for monotonicity re- 
strictions on the intensity function. However, in addition to providing a more general approach 
for intensity estimation, the main feature of this paper is an extension of the intensity mixture 
framework to modeling marked Poisson processes. Indeed, the advantage of a Bayesian non- 
parametric model-based approach will be most clear when it is combined with modeling for 
the conditional mark distribution, thus providing unified inference for point pattern data. 



General theoretical background on Poisson processes can be found, for instance, in Cressie 
( fT993| ), |Klngmanl ( |T993] ), and [Daley and Vere- Jones] ( |2003l ). |Diggle| ( |2003] ) reviews likelihood 



and classical nonparametric inference for spatial NHPPs, and M0ller and Waagepetersen ( 2004 1 
discusses work on simulation-based inference for spatial point processes. 
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A standard approach to (approximate) Bayesian inference for NHPPs is based upon log- 
Gaussian Cox process models, wherein the random intensity function is modeled on logarith- 
mic scale as a Gaussian process (e.g., M0ller et al. 1998[ Brix and Diggle, 200 1[ Brix and| 



M0ller 2001). In particular, Liang et al. (2009) present a Bayesian hierarchical model for 



marked Poisson processes through an extension of the log-Gaussian Cox process to accommo- 
date different types of covariate information. Early Bayesian nonparametric modeling focused 
on the cumulative intensity function, \{s)ds, for temporal point processes, including mod- 
els based on gamma, beta or general Levy process priors (e.g., |Hjort[ |1990[ |Lo| |1992[ |Kuo 



and Ghosh 1997; Gutierrez-Pena and Nieto-Barajas, 2003). An alternative approach is found 



m 



Heikkinen and Arjas (1998 1999), where piecewise constant functions, driven by Voronoi 



tessellations and Markov random field priors, are used to model spatial NHPP intensities. 
The framework considered herein is more closely related to approaches that involve a mix- 



ture model for A(-). In particular, Lo and Weng ( 1989 ) and Ishwaran and James (2004) utilize 
a mixture representation for the intensity function based upon a convolution of non-negative 



kernels with a weighted gamma process. Moreover, Wolpert and Ickstadt ( 1998 ) include the 
gamma process as a special case of convolutions with a general Levy random field, while 



Ickstadt and Wolpert (1999) and Best et al. (2000) describe extensions of the gamma process 



convolution model to regression settings. [Ickstadt and Wolpert| ( |1999| ) also provide a connec- 
tion to modeling for marked processes through an additive intensity formulation. Since these 
mixture models have the integrated intensity term linked to their nonparametric prior for A(-), 
they can be cast as a generalization of our model of independent An. 

A distinguishing feature of the proposed approach is that it builds the modeling from the 
NHPP density. By casting the nonparametric modeling component as a density estimation 
problem, we can develop flexible classes of nonparametric mixture models that allow rela- 
tively easy prior specification and posterior simulation, and enable modeling for multivariate 
mark distributions comprising both categorical and continuous marks. Most importantly, in 
the context of marked NHPPs, the methodology proposed herein provides a unified inference 
framework for the joint location-mark process, the marginal point process, and the conditional 
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mark distribution. In this way, our framework offers a nice simplification of some of the more 
general models discussed in the literature, providing an easily interpretable platform for applied 
inference about marked Poisson processes. The combination of model flexibility and relative 
simplicity of our approach stands in contrast to various extensions of Gaussian process frame- 
works: continuous marks lead to additional correlation function modeling or a separate mark 
distribution model; it is not trivial to incorporate categorical marks; and a spatially changing 
intensity surface requires complicated non-stationary spatial correlation. 

The plan for the paper is as follows. Section |2] presents our general framework of model 
specification for the intensity function of unmarked temporal or spatial NHPPs. Section |3] ex- 
tends the modeling framework to general marked Poisson processes in both a semiparametric 
and fully nonparametric manner. Section |4] contains the necessary details for application of 
the models developed in Sections [2] and [3| including posterior simulation and inference, prior 
specification, and model checking (with some of the technical details given in an Appendix). 



We note that Section 4.2 discusses general methodology related to conditional inference under 
a DP mixture model framework, and is thus relevant beyond the application to NHPP model- 
ing. Finally, Section |5] illustrates the methodology through three data examples, and Section [6] 
concludes with discussion. 



2 Mixture specification for process intensity 

This section outlines the various models for unmarked NHPPs which underlie our general 
framework. As described in the introduction, the ability to factor the likelihood as in ([!]) al- 
lows for modeling of /(x) = A(x)/A7^, the process density, independent of Is.-ji, the integrated 
process intensity. The Poisson assumption implies that is sufficient for At^ in the poste- 
rior distribution and, in Section 4, we describe standard inference under both conjugate and 
reference priors for A?^. Because the process density has domain restricted to the observation 
window IZ, we seek flexible models for densities with bounded support that can provide infer- 
ence for the NHPP intensity and its functionals without relying on specific parametric forms or 
asymptotic arguments. 
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We propose a general family of models for NHPP densities /(x) built through DP mixtures 
/(x; G) of arbitrary kernels, k''(x; 9), with support on TZ. Specifically, 

/(x; G) = J k''(x; e)dG{e), with k''(x; 0) = for x ^ 7^, and G ~ DP(a, Gq), (2) 

where 9 is the (typically multi-dimensional) kernel parametrization. The kernel support re- 
striction guarantees that /(x; G)dx = 1 and hence At^ = A(x)(ix. The random mixing 
distribution G is assigned a DP prior ( Ferguson] [1973 ; Antoniak, 1974[ ) with precision parame- 



ter a and base (centering) distribution Gq{-) = Go{-; ip) which depends on hyperparameters 
For later reference, recall the DP constructive definition ( Sethuraman[ 1994[ ) according to which 



the DP generates (almost surely) discrete distributions with a countable number of atoms drawn 
i.i.d. from Gq. The corresponding weights are generated using a stick-breaking mechanism 
based on i.i.d. Beta(l, a) (a beta distribution with mean (1 + a)^^) draws, {^^ : s = 1, 2, ...} 
(drawn independently of the atoms); specifically, the first weight is equal to Ci and, for / > 2, 
the /-th weight is given by nl=i ~ Cs) ■ The choice of a DP prior allows us to draw from the 
existing theory, and to utilize well-established techniques for simulation-based model fitting. 

The remainder of this section describes options for specification of the kernel and base 
distribution for the model in ([2]): for temporal processes in Section 2. 1 and for spatial processes 
in Section 2.2. In full generality, NHPPs may be defined over an unbounded space, so long as 
the intensity is locally integrable, but in most applications the observation window is bounded 
and this will be a characteristic of our modeling framework. Indeed, the specification of DP 
mixture models for densities with bounded support is a useful aspect of this work in its own 
right. Hence, temporal point processes can be rescaled to the unit interval, and we will thus 
assume that TZ = (0,1). Furthermore, we assume that spatial processes are observed over 
rectangular support, such that the observation window can also be rescaled, in particular, IZ = 



(0, 1) X (0, 1) in Section 2.2 and elsewhere for spatial data. 
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2. 1 Temporal Poisson processes 

Denote by {ti, . . . , tAr} the temporal point pattern observed in interval TZ = (0,1), after the 
rescaling described above. Following our factorization of the intensity as A(t) = A-jif^t) and 
conditional on N, the observations are assumed to arise i.i.d. from f{t; G) = J k^{t; 9)dG{6) 
and G is assigned a DP prior as in We next consider specification for k*(t; 6). 

Noting that mixtures of beta densities can approximate arbitrarily well any continuous den- 
sity defined on a bounded interval (e.g., Diaconis and Ylvisaker[ 1985, Theorem 1), the beta 



emerges as a natural choice for the NHPP density kernel. Therefore, the DP mixture of beta 
densities model for the NHPP intensity is given by 



X{t;G) = An j Kt;fx,T)dG{fx,T), tG(0,l); G~DP(a,Go). 



(3) 



Here, b(-; /i, r) denotes the density of the beta distribution parametrized in terms of its mean 
/X G (0, 1) and a scale parameter r > 0, i.e., b{t; fi, r) oc t^'^-^{l - t)^(i-^)-\ t G (0, 1). Re- 
garding the DP centering distribution Go = Goif^, t), we work with independent components, 
specifically, a uniform distribution on (0, 1) for fi, and an inverse gamma distribution for r with 
fixed shape parameter c and mean (3/{c — 1) (provided c > 1). Hence, the density of Gq is 
5'o(/^, t) oc r"^'^"'"^) exp(— /3r"^)l^g(o,i), where (5 can be assigned an exponential hyperprior. 

The beta kernel is appealing due to its flexibility and the fact that it is directly bounded to 
the unit interval. However, there are no commonly used conjugate priors for its parameters; 
there are conjugate priors for parameters of the exponential family representation of the beta 
density, such as the beta-conjugate distribution in Grunwald et al. (1993]), but none of these are 



easy to work with or intuitive to specify. There are substantial benefits (refer to Section 4) to 
be gained from the Rao-Blackwellization of posterior inference for mixture models (see, e.g.. 



MacEachem et al. 1999 for empirical demonstration of the improvement in estimators) that is 
only possible with conditional conjugacy - that is, in this context, when the base distribution 
is conjugate for the kernel parametrization. Moreover, the nonparametric mixture allows infer- 
ence to be robust to a variety of reasonable kernels, such that the convenience of conjugacy will 
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not usually detract from the quality of analysis. 

We are thus motivated to provide a conditionally conjugate alternative to the beta model, 
and do so by first applying a logit transformation, logit(t) = log — t)), t E (0, 1), and 
then using a Gaussian density kernel. In detail, the logit-normal DP mixture model is then, 

A(t; G) = A7^ y N (logit(t); fi, a^) j^^^dG{fi, a^), te (0, 1); G ~ DP(a, Go). (4) 



The base distribution is taken to be of the standard conjugate form (as in, e.g., Escobar and 
West 1995), such that go{^,(7'^) = N(/i; 5, cr^/fi;)ga(cr"^; z/, w), where ga(-;z/,a;) denotes the 



gamma density with E[(T"^] = v/uj. A gamma prior is placed on u whereas k, v and 5 are fixed 
(however, a normal prior for 5 can be readily added). 

The price paid for conditional conjugacy is that the logit-normal model is susceptible to 
boundary effects: the density specification in (|4]) must be zero in the limit as t approaches the 
boundaries of the observation window (such that logit (t) — > ±oo). In contrast, the beta model 
is not restricted to any single type of boundary behavior, and will thus be more appropriate 
whenever there is a need to model processes which maintain high intensity at the edge of the 
observation window. Section 5 offers empirical comparison of the two models. 

The beta and logit-normal mixtures form the basis for our approach to modeling marked 
Poisson processes, and Section 2.2 will extend these models to spatial NHPPs. Both schemes 
are developed to be as flexible as possible, in accordance with our semiparametric strategy of 
having point event data restricted by the Poisson assumption but modeled with an unrestricted 
NHPP density. However, in some situations it may be of interest to constrain the model further 
by making structural assumptions about the NHPP density, including monotonicity assump- 



tions for the intensity function as in, for example, software reliability applications (e.g., Kuo 



and Yangl 1996). To model monotonic intensities for temporal NHPPs, we can employ the 



representation of non-increasing densities on R+ as scale mixtures of uniform densities. In 
particular, for any non-increasing density h[-) on 1R+ there exists a distribution function G, 



with support on R"*", such that h{t) = h{t\ G) = J 9 ^tte{o,e)<iG{9) (see, e.g., Brunner and Lo 



1989| [Kottas and Gelfandj |2001| ). In the context of NHPPs, a DP mixture formulation could 



be written X{t; G) = An J 0"^lte(o,0)c/G'(^), t G (0, 1), with G ~ DP(a, Go), where Go has 
support on (0, 1), e.g., it can be defined by a beta distribution. Then, A(t; G) defines a prior 
model for non-increasing intensities. Similarly, a prior model for non-decreasing NHPP inten- 
sities can be built from /(t; G) = / 9'^l(^t-i)e(-efl)dG{e), t E (0, 1), with G ~ DP(a, Go), 
where again Go has support on (0, 1). 



2.2 Spatial Poisson processes 



We now present modeling for spatial NHPPs as an extension of the framework in Section 2.1 
As mentioned previously, we assume that the bounded event data has been rescaled such that 
point locations {xi, . . . ,xjv} all lie within the unit square, TZ = (0, 1) x (0, 1). The extra 
implicit assumption of a rectangular observation window is standard in the literature on spatial 



Poisson process modehng (see, e.g., Diggle[ 2003 ). 



The most simple extension of our models for temporal NHPPs is to build a bivariate kernel 
out of two independent densities. For example, a two-dimensional version of the beta mixture 
density in ([3| could be written /(x; G) = / b{xi; fii,Ti)b{x2; fi2,T2)dG{n,T), where n = 
(/ii,/i2) and T = (ti, T2). However, although dependence between xi and X2 will be induced 
by mixing, it will typically be more efficient to allow for explicit dependence in the kernel. A 
possible two-dimensional extension of Q is that of Kottas and Sanso (2007), which employs 



a Sarmanov dependence factor to induce a bounded bivariate density with beta marginals. The 
corresponding model for the spatial NHPP intensity is given by 

A(x;G)=A7^ j b{xi;iii,Ti)b{x2;fi2,r2){l + p{xi-fii){x2-lJ^2))dG{n,T,p), (5) 

where G ~ DP (a, Go) and Go is built from independent centering distributions as in ^ for 
each dimension, multiplied by a conditional uniform distribution for p over the region such 
that 1 + p{xi - pi){x2 - P2) > 0, for all x G 7^. Thus, go{n,T,p) = lpe(c^,C'-)(C"^ - 
C'lj.y^ nLiga(ri"^; z/i,/3,)l^^e(o,i), where G^ = - (max{;Ui/i2, (1 - /ii)(l - 1^2)})'^ and 
= — (min{/ii(/i2 — l),p2iPi — 1)})^- Gamma hyperpriors can be placed on /3i and 132. 
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Model ([5]) has appealing flexibility, including resistance to edge effects, but a lack of condi- 
tional conjugacy requires the use of an augmented Metropolis-Hastings algorithm for posterior 
simulation (discussed in Appendix A.2). The inefficiency of this approach is only confounded 
in higher dimensions, and becomes especially problematic when we extend the models to in- 
corporate process marks. Hence, we are again motivated to seek a conditionally conjugate 
alternative for spatial NHPPs, and this is achieved in a straightforward manner by applying in- 
dividual logit transformations to each coordinate dimension and mixing over bivariate Gaussian 
density kernels. Specifically, the spatial NHPP logit-normal model is 

A(x;G)=A7^ / N(logit(x);^,S)-2 y -t^G(^,S), G~DP(«,G'o), (6) 

where logit(x) is shorthand for [logit(xi), logit(a;2)]'. The base distribution is again of the 
standard conjugate form, such that goifJ', S) = N(^i; S, S/k)W(S^^; u, fl), with fixed n, u, d 
and a Wishart hyperprior for fl. Here, W(-; u, fl) denotes a Wishart density such that E[S^^] = 
ufl'^ andE[S] = - ^y^fl. 

3 Frameworks for modeling marked Poisson processes 

The models for unmarked NHPPs, as introduced in Section 2, are essentially density estimators 
for distributions with bounded support. As mentioned in the Introduction, the Bayesian non- 
parametric approach is most powerful when embedded in a more complex model for marked 
point processes. Section 3.1 describes how the methodology of Section 2 can be coupled with 
general regression modeling for marks, whereas in Section 3.2, we develop a fully nonpara- 
metric Bayesian modeling framework for marked Poisson processes. 

3.1 Semiparametric modeling for the mark distribution 

In the standard marked point process setting, one is interested in inference for the process 
intensity over time or space and the associated conditional distribution for the marks. 

Regarding the data structure, for each temporal or spatial point Xj, i = 1,...,N, in the 
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observation window 71 there is an associated mark taking values in the mark space Ai, 
which may be multivariate and may comprise both categorical and continuous variables. Let 
/i(y I x) denote the conditional mark density at point x. (Note that we use y and y^ as simplified 
notation for y(x) and y(xj).) Under the semiparametric approach, we build the joint model for 
the marks and the point process intensity through 

0(x, y) = A(x)/i(y | x) = Anf{^)h{y | x), x G 7^, y G A^. (7) 



Note that the conditioning in /i(y | x) does not involve any portion of the point process other 
than point x; for instance, for temporal processes, the conditional mark density at time t does 
not depend on earlier times t' < t. Under this setting, the Marking theorem (e.g., proposition 
3.9 in |M0ller and Waagepetersen[ |2004| |Kingmanl |1993[ p. 55) yields that the marked point 



process {(x, y) : x G 7^, y G A^} is a NHPP with intensity function given by (jv]) for (x, y) G 
TZ X A4, and by its extension to i3 x for any bounded B D TZ. 

This intensity factorization, combined with the general NHPP likelihood factorization in 
([1]), results in convenient semiparametric modeling formulations for the marked process through 
a DP mixture model for /(•) (as in Section [2]) and a separate parametric or semiparametric re- 
gression specification for the conditional mark distribution. In particular, assuming that the 
marks {yi}^i are mutually independent given {x^}^]^, and combining ([!]) and ([v]), we obtain 



N N 



P ({x„ y.}^i; An, /(■), H')) oc exp{-An) H fi^i) U Ky. I x.), (8) 



i=l i=l 



such that the conditional mark density can be modeled independent of the process intensity. 

The consequence of this factorization of integrated intensity, process density, and the condi- 
tional mark density, is that any regression model for h can be added onto the modeling schemes 
of Section [2] and provide an extension to marked processes. In some applications, it will be 
desirable to use flexible semiparametric specifications for h, such as a Gaussian process regres- 
sion model, while in other settings it will be useful to fit h parametrically, such as through the 
use of a generalized linear model. As an illustration. Section |5T explores a Gaussian process- 
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based specification, however, the important point is that this aspect of the modeling does not 
require any further development of the underlying nonparametric model for the NHPP inten- 
sity. Moreover, despite the posterior independence of / and h, combining them as in Q leads 
to a practical semiparametric inference framework for the joint mark-location Poisson process. 
The fully nonparametric approach developed in the following section provides an alternative 
for settings where further modeling flexibility is needed. 

3.2 Fully nonparametric joint and implied conditional mark modeling 

While the semiparametric approach of Section 3 . 1 provides a convenient extension of the NHPP 
models in Section 2, the connection between joint and marked processes provides the oppor- 
tunity to build fully nonparametric models for marked point event data. Here, we introduce a 
general modeling approach, built through fully nonparametric models for joint mark-location 
Poisson processes, and describe how this provides a unified inference framework for the joint 
process, the conditional mark distribution, and the marginal point process. 

Instead of specifying directly a model for the marked process, we begin by writing the joint 
Poisson process, PoP(7^ x Ai,(p), defined over the joint location-mark observation window 
with intensity 0(x, y). The inverse of the marking theorem used to obtain equation ^ holds 
that, if the marginal intensity 0(x, y)(iy = A(x) is locally integrable, then the joint process 
just defined is also the marked Poisson process of interest. 

Analogously to the model development in Section 2, we define a process over the joint 
location-mark space with intensity function 

0(x, y; G) = AnJ k^(x; e")ky(y; ey)dG{9^, 9^) = A^fi^, y; G), G ~ DP(a, Go), (9) 

where the mark kernel {y]9^) has support on Ai and the integrated intensity can be defined in 
terms of either the joint or marginal process, such that At^ = J^A(x)c/x = [J^ 0(x, y)c?y] (ix, 
Note that the marginal intensity, and hence the marked point process, are properly defined with 
locally integrable intensity functions. Specifically, we can move integration over Ai inside the 
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infinite sum and 



/ 0(x,y)rfy = An [ k'^lx; ^ / 



UM 



c/G(r,0y) (10) 



An j ]<^^{^;e^)dG^{e^) = Anfi^;G) = A(x). 



Here, G^{9^) is the marginal mixing distribution, which has an implied DP prior with base 
density gQ{9^) = f goiO^, 0^)d6'^ , and we have thus recovered the original DP mixture model 
of Section 2 for the marginal location NHPP PoP(7^, A). As an aside we note that, through 
a similar argument and since 0(x, y) = A(x)/i(y | x), the joint location-mark process of ^ 



satisfies the requirements of proposition 3.9 in |M0ller and Waagepetersen| ( |2004| ), and hence 



the marks alone are marginally distributed as a Poisson process defined on M. with intensity 

0(x, y)rfx = Anj ky(y; ey)dGy{ey). 

In general, both the mixture kernel and base distributions will be built from independent 
components corresponding to marks and to locations, and the random mixing measure is relied 
upon to induce dependence between these random variables. This technique has been employed 



in regression settings by Taddy and Kottas (2010), and provides a fairly automatic procedure 



for nonparametric model building in mixed data-type settings. For example, suppose that a 
spatial point process is accompanied by categorical marks, such that marks . . . , un} are 
each a member of the set = {1,2,..., M}. The joint intensity model can be specified as 

<t>{^,y-G) = An j k'^(x; q), G ~ DP(«, G-(r)Dir(q; a)), (11) 

where q = [gi, . . . , ^m] is a probability vector with qy = Ft(Y = y \ q), Dir(q; a) is the 
Dirichlet distribution, with a = (ai, um), such that E,{qy \ a) = cby/ Y^^iO-s, and the 
location-specific kernel, k^, and centering distribution, Gq, are specified as in either (|5]) or 
(|6]) and thereafter. Additional marks can be incorporated in the same manner by including 
additional independent kernel and base distribution components. 

Similarly, continuous marks can be modeled through an appropriate choice for the indepen- 
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dent mark kernel. For example, in the case of real-valued continuous marks (i.e., = M) for 
a temporal point process, the choice of a normal density kernel leads to the intensity model 

y- G) = AnJ k\t; r/, a')dGie\ r/, a'), DP («, a^)) . (12) 



The location specific kernel, k*, and base measure, G\, can be taken from Section 2.1 , Gq can 
be specified through the conjugate normal inverse-gamma form as in ([4]). Other possible mark 
kernels are negative-binomial or Poisson for count data (as in Section [S!2] ), a WeibuU for failure 



time data, or a log-normal for positive continuous marks (as in Section 5.3 ). 

As an alternative to this generic independent kernel approach, the special case of a com- 
bination of real-valued continuous marks with the logit-normal kernel models in either (|4]) or 



(|6]) allows for joint multivariate-normal kernels. Thus, instead of the model in ( 12), a temporal 



point process with continuous marks is specified via bivariate normal kernels as 

0(t,y;G) = A7^ j N ([logit(t), y]'; ^, S) S), G~DP(«,G'o), (13) 

with base distribution of the standard conjugate form, exactly as described following (|6]). Spec- 
ification is easily adapted to spatial processes or multivariate continuous marks through the use 



of higher dimensional normal kernels (see Section 5.3 for an illustration). 

A key feature of the joint mixture modeling framework for the location-mark process is 
that it can provide flexible specifications for multivariate mark distributions comprising both 
categorical and continuous marks. For any of the joint intensity models specified in this section, 
inference for the conditional mark density is available through 

ftiy|x,aj- ^^^.^^ - jk-(x;^^-)dG'-(e-) • ^ ^ 

Of course, other conditioning arguments are also possible if, for example, some subset of the 



marks is viewed as covariates for a specific mark of interest. In any case, the integrals in ( 14) 



are actually infinite sums induced by discrete realizations from the posterior distribution for G. 
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In Section 4.2. we show that truncation approximations to the infinite sums allow for proper 
conditional inference and, hence, for fully nonparametric inference about any functional of the 
conditional mark distribution. 



4 Implementation 

This section provides guidelines for application of the models proposed in Sections |2] and |3} 



with prior specification and posterior simulation briefly discussed in Section 4.1 (further details 



can be found in the Appendix), inference for marked NHPP functionals in Section 4.2, and 



model checking in Section 4.3 



4. 1 Prior specification and posterior simulation 

As with our approach to model building, we can specify the prior for integrated intensity inde- 
pendent of the prior for parameters of the DP mixture density model. The marginal likelihood 
for At^ corresponds to a Poisson density for A^, such that the conjugate prior for K-ji is a gamma 
distribution. As a default alternative, we make use of the (improper) reference prior for At^, 
which can be derived as Ti{!s.n) oc A^^ for At^ > 0. The posterior distribution for the integrated 
intensity is then available analytically as a gamma distribution, since the posterior distribution 
for the NHPP intensity factorizes as p (/(•), At^ | data) = p(/(-) | data)p(A7^ | A^). In 
particular, p(A7^ | A^) = ga(A^, 1) under our default reference prior. Similarly, under the semi- 



parametric approach of Section 3.1 prior specification and posterior inference for any model 
applied to the conditional mark distribution can be dealt with separately from the intensity 
function model, and will generally draw on existing techniques for the regression model of 
interest. 

What remains is to establish general prior specification and MCMC simulation algorithms 
for the DP mixture process density models of Sections |2]and 3.2 In a major benefit of our 



approach - one which should facilitate application of these models - we are able here to make 
use of standard results and methodology from the large literature on DP mixture models. Our 
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practical implementation guidelines are detailed in the Appendix, with prior specification in 
A.l and a posterior simulation framework in A. 2. 



4.2 Inference about NHPP junctionals 

Here, we describe the methods for posterior inference about joint or marginal intensity func- 
tions and for conditional density functions. We outline inference for a general NHPP with 
events {zj}^^, possibly consisting of both point location and marks, and leave specifics to the 
examples of Section |5j 

Due to the almost sure discreteness of the DP, a generic representation for the various mix- 
ture models for NHPP densities is given by /(z; G) = J2'^^pi'k(z;'di), where the 'di, given 
the base distribution hyperparameters 'ip, are i.i.d. from Gq, and the weights pi are generated 
according to the stick-breaking process discussed in Section [2} Here, z may include only point 
locations (as in the models of Section [5]) or both point locations and marks whence k(z; = 



k''(x; t9^)ky(y; ^9^) (as in Section 3.2). Hence, the DP induces a clustering of observations: 



for data = {zi,...,Ziv}, ifwe introduce latent mixing parameters = {9i, . . . ,9^} such that 
z, I Oi k(z,; Oi), with 0i I G ~ G, for i = 1, . . . , iV, and G I a, ~ DP(a, Go{-, ip)), then 
observations can be grouped according to the number, m < A^, of distinct mixing parameters 
in 0. This group of distinct parameter sets, 0* = {9^, . . . , 9^}, maps back to data through 
the latent allocation vector, s = [si, . . . , sat], such that 9i = 9*.. The expanded parametriza- 
tion is completed by the number of observations allocated to each unique component, n = 
[ni, . . . , rijn], where Uj = ^s,=j, and the associated groups of observations {zj : Si = j}. 
If G is marginalized over its DP prior, we obtain the Polya urn expression for the DP prior 
predictive distribution, 

m 

p{9o I 0\a,^) = dE[G{9o) \ 0\a,^] (x ago{9o;^) + Y,^,Se*i^o) (15) 

where 5a denotes a point mass at a. Moreover, based on the DP Polya urn structure, the prior 
for 0*, given m and ip, is such that 9* \ ip Go{-', V^), for j = 1, . . • , m. 
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Within the DP mixture framework, estimation of linear functionals of the mixture is possible 
via posterior expectations conditional on only this finite dimensional representation (i.e., it is 
not necessary to draw G). In particular, with the NHPP density modeled as our generic DP mix- 
ture, the posterior expectation for the intensity function can be written as E [A(z; G) | data] = 
E(A7^ I A^)p(z I data), where p(z | data) = E [/(z; G) | data] is the posterior predictive density 



given by 



a + N 



a 



k(z; e)dGo{9] ip) + njk{z; 9*) ) p(r , s,a,ip \ data)de*dsdadip. (16) 



Hence, a point estimate for the intensity function is readily available through E [/(z; G) \ data] 
estimated as the average, for each point in a grid in z, over realizations of ([16]) calculated for 
each MCMC posterior sample for 0*, s, a and ip- 

However, care must be taken when moving to posterior inference about the conditional 



mark distribution in (14). As a general point on conditioning in DP mixture models for joint 



distributions, Polya um-based posterior expectation calculations, such as (16), are invalid for 



the estimation of non-linear functionals of A or /. For example, Miiller et al. ( 1996) develop 
a DP mixture curve fitting approach that, in the context of our model, would estimate the 
conditional mark density by 

= J /kx(x;gx),E[G(g)|^,a,^] P^^'"'^ " d.t.)dedadi., (17) 
which is the ratio of Polya urn joint and marginal density point estimates given 6 and DP prior 



parameters a, ip, averaged over MCMC draws for these parameters. Unfortunately, ( |T7[ ) is not 
E [h{y I x; G) | data], the posterior expectation for random conditional density h{y | x; G) = 
/(x, y; G)/ fix; G), which would be the natural estimate for the conditional mark density at 



any specified combination of values (x, y). Hence, the regression estimate in Miiller et al 



( |1996 ) as well as that proposed in the more recent work of Rodriguez et al. ( 2009| ), based on 



p(x, y I data)/p(x | data), provide only approximations to E [/i(y | x; G) | data]; in particular, 
the latter estimate is approximating the expectation of a ratio with the ratio of expectations. 
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Such approximations become particularly difficult to justify if one seeks inference for non- 
linear functional of /i(y | x; G). 

Hence, to obtain the exact point estimate E [h{'y | x; G) | data], and, most importantly, to 
quantify full posterior uncertainty about general functional of the NHPP intensity, it is neces- 
sary to obtain posterior samples for the mixing distribution, G. Note that p(G | data) = J p(G | 
0*^, s, a, 4))'p{0*, s, a, ^ | data) d0*dsda dip, where p(G | 0*, s, a, ip) follows a DP distribution 
with precision parameter a + N and base distribution given by ( [T5| ) (see Appendix A. 2). As 
discussed in Ishwaran and Zarepour ( |2002[ ), using results from Pitman ( 1996 1, a draw for G \ 
6* , s, a, can be represented as goG^*(") + YlJ=i Qj^e^i-), where G* \ a,ip ^ DP(q;, Go{'ijj)), 
and, independently of G*, (go, gi, g™) | a,s ~ Dir(go, gi, gm; «, ^m)- Therefore, 

posterior realizations for G can be efficiently generated, by drawing for each posterior sample 
{r,s,a,^}, 

{L "I m 

1=1 J J=l 

that is, using a truncation approximation to G* based on the DP stick-breaking definition. 
Specifically, the i^i, I = 1, ...,L, are i.i.d. from G'o(^/'), and the pi are constructed through 
i.i.d. Beta(l,a) draws, (s, s = 1,...,L - 1, such that pi = d, pi = Oni=i(l " Cs), for 
I = 2, L — 1, and p^ = 1 ~ J2f=i Pi- The truncation level L can be chosen using standard 
distributional properties for the weights in the DP representation for G* = Xlj^i For 
instance, E(^^^ | a) = 1 — {a/(« + l)}^, which can be averaged over the prior for « to es- 
timate IE(X]^i '^i)- Given any specified tolerance level for the approximation, this expression 
yields the corresponding value L. Note that even for dispersed priors for a, relatively small 
values for L (i.e., around 50) will generally provide very accurate truncation approximations. 

Now, the posterior distribution for any functional (linear or non-linear) of the NHPP den- 
sity, and thus of the intensity function, can be sampled by evaluating the functional using the 
posterior realizations Gl- For example, suppose that z = [t, y], such that we have a temporal 
process with a single mark, where the mixture kernel factors as k(z; 9) = k^{t; 6'*)k^(?/; 9^). 
Given a posterior realization for Gl and a posterior draw for A-ji, a posterior realization for the 



18 



marginal process intensity at time t is available as 



A(t; Gl) = An [go Eti^'^'^*^' + 5Zj= i ^, 



where 'di = {'d\^'d\) and 0* = (6'**, 6**^), and a realization for the conditional density of mark 
value y at time t arises through 

/ifu t: Gl) = f — . (18) 

Similarly, realized conditional expectation is available as 

Hy I t; Gl] = {fit; GlT' {go Eti^'^*^^' ^'^^^^ ' + 5^ -li ^^'^^^^^ ' ^^''^ 

(19) 

a weighted average of kernel means with time-dependent weights. For multivariate Gaussian 



kernels, as in ( 13), one would use conditional kernel means (available through standard multi- 



variate normal theory; see Section 5.2). The approach applies similarly to multivariate marks 
and/or to marked spatial NHPP, and we can thus obtain flexible inference for general function- 
als of marked NHPPs with full uncertainty quantification. 

4.3 Model checking 

A basic assumption implied by the Poisson process model is that the number of events within 
any subregion of the observation window are Poisson distributed, with mean equal to the in- 
tegrated intensity over that subregion. Hence, a standard approach to assessing model validity 
is to compare observed counts to integrated intensity within a set of (possibly overlapping) 



subregions (e.g., Digglel 2003; Baddeley et al. . 2005). 



An alternative approach to model checking is to look at goodness-of-fit for simplifying 
transformations of the observations. In particular, we propose transforming each margin of the 
point event data (i.e., each spatial coordinate and each mark) into quantities that are assumed, 
conditional on the intensity model, distributed as i.i.d. uniform random variables. Posterior 
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samples of these (assumed) i.i.d. uniform sets can be compared, either graphically or formally, 
to the uniform distribution to provide a measure of model validity. 

Consider first temporal point processes, and assume that the point pattern {ti : i = 1,...,A^}, 
with ordered time points = to < ^ h < ••• < ^Af < 1, is a realization from a NHPP with 
intensity function \{t) and cumulative intensity function A{t) = X{s)ds. Then, based on 
the time-rescaling theorem (e.g., Daley and Vere-Jones[ 2003 [ ), the transformed point pattern 
{A(tj) : i = 1, A^} is a realization from a homogeneous Poisson process with unit rate. Let 
A(t; Gl) be the posterior draws for the cumulative intensity, obtained following the approach of 



Section 4.2 Then, with A(0; Gl) = by definition, the rescaled times A(tj; Gl) — A(tj_i; Gl), 
i = 1, N, are independent exponential random variables with mean one. Thus, the sampled 
Ui = 1 — exp{ — (A(tj; Gl) — A(tj_i; Gl))}, i = 1, A^, are independent uniform random 
variables on (0, 1). 

This approach can be extended to spatial processes by applying the rescaling to each mar- 
gin of the observation window (e.g., Cressie[ 1993 1. If we have data corresponding to a NHPP 
on 7^ = (0, 1) X (0, 1) with intensity A(x), then point event locations along (say) the first 
margin of the window are the realization of a one-dimensional NHPP with intensity Ai(xi) = 

X{x.)dx2, and analogously for A2(x2). Since the kernels in (jsj) and (|6) are easily marginal- 
ized, cumulative intensities Ai(-) and A2(-) are straightforward to calculate as sums of marginal 
kernel distribution functions, based on the sampled Gl as described in Section 4.2. For each 
dimension j, these are then applied to ordered marginals {xj i, . . . , xj^n} to obtain i.i.d. uni- 
form random variables, Uij = 1 — exp{ — (Aj(xj,j; Gl) — Aj{xj-i^i; Gl))}, « = 1, A^, where 
by definition Aj{xjfi; Gl) = for j = 1, 2. 

Finally, there are a variety of ways that the marks can be transformed into uniform random 
variables (for instance, the marginal process for continuous marks is also Poisson, such that the 
time-rescaling theorem applies), but, arguably, the most informative approach is to look at the 



conditional mark distribution of ( fT4| ). Full inference is available for the conditional cumulative 
distribution function H(y \ x; G^) = J^^h(s \ x;GL)ds, through a summation similar to 



that in ( 18 1, at any desired points (x, y). We thus obtain sets of Ui that are assumed to be i.i.d. 
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uniform by taking, for each sampled Gl, the distribution function evaluated at the data such 
that Ui = H{yi \ Xj; Gl), for i = 1, . . . , A^. 

Goodness-of-fit is evaluated through comparison of the Ui samples with the uniform dis- 
tribution, using either graphical or distance-based techniques. For instance, in the context of 



neuronal data analysis. Brown et al. (2001) used standard tests and quantile-quantile (Q-Q) 
plots to measure agreement of the estimated Ui with the uniform distribution on (0, 1). In the 
examples of Section [5| we focus on Q-Q plots for graphical model assessment, and find that 
these provide an intuitive picture of the marginal fit. In particular, under our Bayesian modeling 
approach, inference about model validity can be based on samples from the full posterior for 
each set of Ui, with each realization corresponding to a single draw for Gl, through plots of 
posterior means and uncertainty bounds for the Q-Q graphs. 

The rescaling diagnostics involve a checking of the fit provided by the DP mixture model as 
well as of the Poisson process model assumption, and thus characterize a general nonparamet- 
ric model assessment technique. Note that, in evaluating the model for event-location intensity, 
it is not, in general, feasible under this approach to distinguish the role of the Poisson assump- 
tion from the form of the nonparametric model for the NHPP density. The flexibility of the 
DP mixture modeling framework is useful in this regard, since by allowing general intensity 
shapes to be uncovered by the data, it enables focusing the goodness-of-fit evaluation on the 
NHPP assumption for the point process. Furthermore, all of these goodness-of-fit assessments 
are focused on model validity with respect to marginal processes (although, of course, these are 
implied marginals from a multidimensional fit). It is possible to extend the rescaling approach 
to higher dimensions, by defining a distance metric in the higher dimensional space and evalu- 



ating cumulative intensity functions with respect to this metric (e.g . , |Diggle[ 1 1 990) . However, 
such procedures are considerably more difficult to implement and will need to be designed 
specifically for the application of interest. 
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5 Examples 



We include three data examples to illustrate the methodology. Specifically, Section 5.1 



m- 



volves a simulated data set from a one-dimensional Poisson process with both categorical and 



continuous marks. In Sections 5.2 and 5.3 we consider real data on coal mining disaster events 



occurring in time with count marks, and on spatial tree locations data with trunk-diameter 
marks, respectively. 



5. 1 Simulated events with continuous and binary marks 

We first consider a simulated data set from a temporal Poisson process with observation window 
7^ = (0, 1) and intensity A(t) = 250 1/11, 11) + h{t; 4/7, 7)), such that Kti = 500. The 
simulated point pattern comprises = 481 point events, which are accompanied by binary 
marks z and continuous marks y generated from a joint conditional density h(y, z \ t) = 
h{y I z,t)Vi{z I t). Here, Vi^z = 1 \ t) = and the conditional distribution for y, given 
z and t, is built from y = -10(1 - + e, with e ~ N(0, 1) if z = 0, and e ~ ga(4, 1) if 
z = 1. Hence, the marginal regression function for y given t is non-linear with non-constant 
error variance, and Pr(z = l\t) increases from to 1 over IZ. 

We consider a fully nonparametric DP mixture model consisting of the beta kernel in (|3]) 
for point locations combined with a normal kernel for y and a Bernoulli kernel for z. Hence, 
the full model for the NHPP density is given by 

f{t,y,z-G) = j b{t; fi,TMy;r,,(j))q'{l - qy-'dG{^I,T,r],<P,q), G^DF{a,Go) 

where go{^, r, r], 0, q) = l^e(o,i)ga(r~i; 2, /3^)N(r/; 0, 2O0)ga(0-^; 2, P^)b{q; 0.5, 1). We use 
the reference prior for At^, and for the DP hyperpriors take a ~ ga(2, 1), (3r ~ ga(l, 1/20) 
and ~ ga(l, 1); note that j3r and are the means for r and (f), respectively, under Go- The 
hyperpriors are specified following the guidelines of Appendix A.l, and posterior simulation 
proceeds as outlined in Appendix A. 2. Since the beta kernel specification is non-conjugate, 
we jointly sample parameters and allocation variables with Metropolis-Hasting draws for each 
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(Hi, Ti) and Si given s*^^'^ and (/x*, t*)*^~*\ as in algorithm 5 of Neal (2000). 

Results are shown in Figure[T} In the top row, we see that our methods are able to capture the 
marginal point intensity and general conditional behavior for y and z\ note that the uncertainty 
bounds are based on a full assessment of posterior uncertainty that is made possible through 
use of the truncated approximations to random mixing measure G (as developed in Section 



4.2). We also fit a Gaussian process (GP) regression model to the (t, y) data pairs (using the 



tgp package for R under default parametrization) and, in contrast to our approach based on 



draws from h{y \ t; Gl) as in ( 18 ), the top middle panel shows the GP model's global variance 
as unable to adapt to a wider skewed error distribution for larger t values. 

The middle row of Figure [T] illustrates behavior for a slice of the conditional mark density 
for y,att = 1/2, both marginally and given 2 = or 1. The marginal (left-most) plot shows 
that our model is able to reproduce the skewed response distribution, while the other two plots 
capture conditional response behavior given each value for z. As one would expect, posterior 
uncertainty around the conditional mark density estimates is highest at the transition from nor- 
mal to gamma errors. Finally, posterior inference for model characteristics is illustrated in the 
bottom row of Figure [T| Peaked posteriors for (3^. and show that it is possible to learn about 
hyperparameters of the DP base distribution for both t and y kernel parameters, despite the 
flexibility of a DP mixture. Moreover, based on the posterior distribution for m, we note that 
the near to 500 observations have been shrunk to (on average) 12 distinct mixture components. 



5.2 Temporal Poisson process with count marks 

Our second example involves a standard data set from the literature, the "coal-mining disasters" 
data (e.g., [Andrews and Herzbergl|1985[ p. 53-56). The point pattern is defined by the times (in 
days) of 191 explosions of fire-damp or coal-dust in mines leading to accidents, involving 10 
or more men killed, over a total time period of 40,550 days, from 15 March 1851 to 22 March 
1962. The data marks y are the number of deaths associated with each accident. 

This example will compare two different mixture models for marginal location intensity: 
a "direct" model with beta-Poisson kernels, and a "transformed" model with data mapped to 
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Figure 1: Simulation study results. On top, from left to right, we have posterior mean and 90% interval 
for the marginal intensity A(i; G) (with the true intensity denoted by the grey hne), the data (dark grey for 
z = 1), and posterior 90% predictive intervals based on both h(y \ t; G) (solid lines) and GP regression 
(dotted lines), and posterior mean and 90% intervals for Pr (z = 1 | t ; G) (with the true function denoted 
by the grey line). The middle row has mean and 90% intervals for conditional densities for y at t = 1/2, 
marginalized over z (left panel) and conditional on z (middle and right panels), with true densities 
plotted in grey. Lastly, the bottom row shows posterior samples for fir and (dark grey, with priors in 
the background) and for the number of latent mixture components. 
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and fit via multivariate normal kernels. The first scheme models data directly on its origi- 
nal scale, but requires Metropolis-Hastings augmented MCMC for the beta kernel parameters, 
and dependence between t and y is induced only through G. The second model affords the 
convenience of the collapsed Gibbs sampler and correlated kernels, but on a transformed scale. 

Following our general modeling approach, both models use the reference prior for A-ji and 
assume NHPP density form f{t, y;G) = J k(t, y; 9)dG{9) with G ~ DP(a, Gq) and 7i{a) = 
ga(2, 1). The distinction between the two models is thus limited to choice of kernel and base 
distribution. For the direct model, 

k(t, /i, r, 0) = 6(t;/i,r)Po>io(i/;0), (20) 
5(o(/i,r,0) = 1^6(o,i)ga(r"^;2,/3^)ga(0; 1,1/60), 

where Po>io(?/; (p) is a Poisson density truncated at y = 10, and with n((3r) = ga(l, 1/63). 
This leads to prior expectations E[0] = 60 and E[r] = E = 63 for mean location kernel 
precision (1 + r)/ (yu(l — fi)) ^ 4(1 + 63), which translates to a standard deviation of 1/16. For 
the transformed model, we take y = y — 9.5 and 

. N([log.tW.log(M';M.E) 

yt{l-t) 

go(M,S) = N(/x;(0,2.5)',10S)iy(S-^3,ri), 

with Tr{n) = W(3, diag[10, 20]) for E(S) = 2/3E(fi) = diag[l/5, 1/10] (logit(t) and log(y) 
range in (-5,5) and (-1,6), respectively). Both models were found to be robust to changes in this 
parametrization (e.g., E[0] G [10, 100] and diagonal elements of E[S] in [0.1, 1]). 

Results under both models are shown in Figure [2} In the top left panel, we see that marginal 
process density estimates derived from each model are generally similar, with the normal model 
perhaps more sensitive to data peaks and troughs. There is no noticeable edge effect for either 
model. The Q-Q plot in the bottom left panel shows roughly similar fit with the normal model 
performing slightly better. The top and bottom right panels report inference for the count mark 
conditional mean and distribution Q-Q plot. For the beta-Poisson model, posterior realizations 
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Figure 2: Coal-mining disasters. Mean and 90% intervals for (clockwise from top-left): marginal density 
f{t]Gi,) (with data histogram); conditional expected count E(y | t;GL) (data counts in grey); and 
posterior Q-Q plots for Pr(y < yi \ U^Gl) and A{ti;GL), respectively. 



for E(?/ I t; Gl) are obtained using ( 19 1. The conditional mean calculation for the normal model 
must account for the correlated kernels (and the transformation to y), such that E(?/ | t; Gl) is 



9 
2 



1=1 j=i 



where E[?/ | t,6] = exp [fiy + pcrt^{t — /i^) + 0.5((t^ — p^cr^"^)] with ^ = {fit, fiy) and S 
partitioned into variances (cr^ , cr^) and correlation p. Similarly, uniform quantiles for the con- 
ditional mark distribution under the beta-Poisson model are available as weighted sums of 
Poisson distribution functions, while the normal model calculation for Pr(2/ < yi \ tf, Gl) is as 
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above for E{y \ t;GL), but with E[y \ t,9] replaced by Pr(|/ < yi \ ti;9) = 
^ (^[jji — fly + pa^^{ti — fit)'j (cr^ — p^(T(^^)^^/^) . The estimated conditional mean functions 
are qualitatively different, with the poisson model missing the peak at WWl . Indeed, the corre- 
sponding QQ plot shows that the normal model provides a better fit to this data; we hypothesize 
that this is due to the equality of mean and variance assumed in Poisson kernels, and may be 
fixed by using instead, say, truncated negative binomials. 



5.3 Spatial Poisson process with continuous marks 

Our final example considers the locations and diameters of 584 Longleaf pine trees in a 200 x 
200 meter patch of forest in Thomas County, GA. The trees were surveyed in 1979 and the 
measured mark is diameter at breast height (1.5 m), or dbh, recorded only for trees with greater 
than 2 cm dbh. The data, available as part of the spats tat package for R, were analyzed 



by Rathbum and Cressie (1994) as part of a space-time survival point process. Poisson pro- 
cesses are generally viewed as an inadequate model for forest patterns, due to the dependent 
birth process by which trees occur. However, the NHPP should be flexible enough to account 
for variability in tree counts at a single time point and, in this example, we will concentrate 
primarily on inference for the conditional dbh mark distribution. 



To analyse this data set, we employ a spatial version of the model in ( 13), with tree marks 



log-transformed to lie on the real line. Thus, our three-dimensional normal kernel model is 
G) ^ A, / N(|log.t(x)^ogfa 2)|-;^.E) ^^^^^ ^ ^ ^^^^^ 

>^ [y [ U=i ^'i) 

The base distribution is taken to be fi-ol/^, S) = N(^i; (0, 0, 1)', 100S)W(S~^; 4, 17), with 
7r(f2) = W(4, diag[0.1, 0.1, 0.1, 0.1]). A ga(2, 1) prior is placed on a. Posterior sampling 
follows the fully collapsed Gibbs algorithm of Appendix A. 2. 

In this data set, high density clusters of juveniles trees {dbh < hem) combine with the more 
even dispersal of larger trees to form conditional mark densities with non-standard shapes and 
non-homogeneous variability. This behavior is clearly exhibited in the posterior estimates of 
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Figure 3: Longleaf pines. The left panel has data (point size proportional to tree diameter) and a Q-Q 
plot (mean and 90% interval) for P h{s \ x; GL)ds evaluated at data. The right panel plots posterior 
mean and 90% intervals for h{y \ x; Gl) at four specific x values. 

the conditional density for dhh, shown on the right side of Figure [3} at four different locations 
in the observations window. Although conditional densities vary in shape over the different 
locations, each appears to show the mixture of a diffuse component for mature trees combined 
with a sharp increase in density at low dhh values, corresponding to collections of juvenile trees 
(only some of whom make it to maturity). It is notable that we are able to infer this structure 
nonparametrically, in contrast to existing approaches where the effect of a tree-age threshold is 
assumed a priori (as in Rathburn and Cressie} 1994). Finally, the conditional mark distribution 
Q-Q plot on the bottom right panel of Figure[3](based on calculations similar to those in Section 
5.2| ) shows a generally decent mean-fit with wide uncertainty bands corresponding to the 95% 
and 5% density percentile Q-Q plots. 



6 Discussion 

We have presented a nonparametric Bayesian modeling framework for marked non-homogeneous 
Poisson processes. The key feature of the approach is that it develops the modeling from the 
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Poisson process density. We have considered various forms of Dirichlet process mixture mod- 
els for this density which, when extended to the joint mark-location process, result in highly 
flexible nonparametric inference for the location intensity as well as for the conditional mark 
distribution. The approach enables modeling and inference for multivariate mark distributions 
comprising both categorical and continuous marks, and is especially appealing with regard to 
the relative simplicity with which it can accommodate spatially correlated marks. We have dis- 
cussed methods for prior specification, posterior simulation and inference, and model checking. 
Finally, three data examples were used to illustrate the proposed methodology. 

The Poisson assumption for marked point processes is what enables us to separate mod- 
ehng for the process density from the integrated intensity. This simplification is particularly 
useful for applications involving several related intensity functions and mark distributions, and 
is less restrictive than it may at first appear. For instance, Taddy ( 2010[ ) presents an estimation 
of weekly violent crime intensity surfaces, using autoregressive modeling for marked spatial 
NHPPs, and |Kottas et al. (2011 ) compares neuronal firing intensities recorded under multiple 
experimental conditions, using hierarchically dependent modeling for temporal NHPPs. 

Among the possible ways to relax the restrictions of the Poisson assumption, while retain- 
ing the appealing structure of the NHPP likelihood, we note the class of multiplicative intensity 
models studied, for instance, in Ishwaran and James ( 2004[ ). These models for marked point 
processes are under the NHPP setting and, indeed, follow the simpler strategy of separate mod- 
ehng for the process intensity and mark density as in the semiparametric framework of Section 
3.1. More generally, one could envision relaxing the Poisson assumption for the number of 
marks through a joint intensity function such that the location intensity is not the marginal of 
the joint intensity over marks. Such extensions would however sacrifice the main feature of our 
proposed framework - flexible modeling for multivariate mark distributions under a practical 
posterior simulation inference scheme. As a more basic extension, our factorization in (1) could 
be combined with alternative specifications for integrated intensity; for example, hierarchical 
models may be useful to connect intensity across observation windows. 
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Appendix: Implementation Details for Dirichlet Process Mixture Models 

A. 1 Prior Specification 

Prior specification for the DP precision parameter is facilitated by the role of a in controlling 
the number, m < A^, of distinct mixture components. For instance, for moderately large A^, 
E[m I a] ^ a log {{a + N)/a). Furthermore, it is common to assume a gamma prior for a, 
such that n{a) = ga{a; a^, ba), and use prior intuition about m combined with E[m | a] to 
guide the choice of aa and ba- 

Specification of the base distribution parameters will clearly depend on kernel choice and 
application details, and DP mixture models are typically robust to reasonable changes in this 
specification. First, the base distribution for kernel location (usually the mean, but possibly 
median) can be specified through a prior guess for the data center; for example, this value 
can be used to fix the mean parameter 5 in ([6]) or the mean of a normal hyperprior for 6. In 
choosing dispersion parameters, note that the DP prior will place most mass on a small number 
of mixture components, with the remaining components assigned very little weight and, hence, 
very few observations. At the same time, this behavior can be overcome in the posterior and 
it is important to not restrict the mixture to overly-dispersed kernels. Thus, the expectation of 
the kernel variance (or scale, or shape) parameters should be specified with a small number of 
mixture components in mind, but with low precision. For example, again in the context of ([6]), 
the square -root of the hyperprior expectation for diagonal elements of fl can be set at 1/8 to 
1/16 of a prior guess at data range, and the precision v will be as small as is practical (usually 
the dimension of the kernel plus 2). The factor k is then chosen to scale the mixture to expected 
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dispersion in fi. 

Moreover, except when specific prior information about co-dependence is available, it is 
best to center Go on kernel parametrization that implies independence between variables, such 
that the mixture is centered on a model with dependence induced nonparametrically by G. 
For example, in the model of ([6]), we assume zeros in the off-diagonal elements for the prior 
expectation of fl, and this is combined with a small u to allow for within-kernel dependence 
where appropriate. A prior expectation of independence also fits with our general approach of 
building kernels for mixed-type data as the product of multiple independent densities. 

Note that we have chosen to introduce prior information into the base measure based on the 
intuition arising from a small number of large mixture components and a near zero. Recent 
work in Bush, Lee, and MacEachem ( 2010| ) provides a rigorous treatment of non-informative 



prior specification, and they advocate a hierarchical scheme for a\Go that maintains desirable 
properties at all scales of precision. As the main work here - use of mixtures for modeling joint 
location-mark Poisson process densities - is independent of prior and base measure choice, 
these innovations, as well as application-specific prior schemes, could potentially be integrated 
into our framework. 

A.2 Posterior Simulation 



Using results from Antoniak (1974]), the posterior distribution for the DP mixture model is 



partitioned as p(G', 0*, s, a, | data) = p(G | 6* , s, a, ?/')p(0*, s^a^ip \ data), where G, given 
6* , s, a, iIj, is distributed as a DP with precision parameter a + N and base distribution given 
by ([15]). Hence, full posterior inference involves sampling for the finite dimensional portion 
of the parameter vector, which is next supplemented with draws from the conditional posterior 
distribution for G (obtained as discussed in Section]?^). A generic Gibbs sampler for posterior 



simulation from p(0*, s, a, ^ | data), derived by combining MCMC methods from MacEachern 
( |1994[ ) and [Escobar and West| ( |1995p , proceeds iteratively as follows: 



For i = 1, . . . ,N , denote by s*- *^ the allocation vector with component removed, and 
by A^i *^ the number of elements of s*^^*) that are equal to s. Then, if s = s,- for some 
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r 7^ the i-th allocation variable is updated according to 

Pr(si = s I s(-*),a,V',data) oc — ^ / kfz^; rwr I s^"*), ^, data)rfr, 

N - 1 + a J 

where p{9* \ s'^~'^\^p, data) is the density proportional to go{9*:, ip) Y[{rj^rsr=s} k(zr; d*)- 
Moreover, the probability of generating a new component, that is, Pr(sj 7^ Sr for allr 7^ 
i I s^~*\ a, ^z^, data), is proportional to a{N — 1 + a)^^ J k(zj; 6*)go{6*; ilj)d6*. 

For j = 1, m, draw 9* from | s, ijj, data) oc 5(0(6*^*; ^) ll{i:s,=j} H^i'^ ^j)- 

Draw the base distribution hyperparameters from 7r(?/;) YYjLi 5'o(^j ; V^)^ where 7r(?/') is the 
prior for 'ip. Finally, if a is assigned a gamma hyperprior, it can be updated conditional 



on only m and N using the auxiliary variable method from Escobar and West ( 1995 1. 



The integrals that are needed to update the components of s can be evaluated analytically for 
models where Go is conjugate for k(-;9). It is for this reason that conditionally conjugate 
mixture models can lead to substantially more efficient posterior sampling, especially when 9 
is high-dimensional. When this is not true (as for, e.g., beta kernel models or the truncated 



Poisson of Equation 20), the draw for s requires use of the auxiliary parameters, 6*, sampled 
as in the second step of our algorithm, in conjunction with a joint Metropolis-Hastings draw 
for each 9i and Sj given 6^^^^ and s*^^*). In particular, we can make use of algorithms from 



Neal 



(2000) for non-conjugate models. 
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